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Abstract. - The response of a nonlinear stochastic system driven by an external sinusoidal 
time dependent force is studied by a variety of numerical and analytical approximations. The 
validity of linear response theory is put to a critical test by comparing its predictions with 
numerical solutions over an extended parameter regime of driving amplitudes and frequencies. 
The relevance of the driving frequency for the applicability of linear response theory is explored. 



The response of dissipative physical systems to small amplitude external perturbations is 
usually described with the powerful tools of linear response theory (LRT) jlj, as it is generally 
accepted that the effect of the perturbation can be described in terms of small deviations from 
the behavior of the unperturbed system. In particular, for long times and for systems which 
in the absence of driving reach an equilibrium distribution, LRT provides an approximate 
expression for the probability distribution obtained by keeping just the linear terms in a series 
expansion in the external amplitude. The purpose of the present letter is to point out the 
relevance of parameters other than the amplitude of the driving force, for the validity of LRT. 
We will show that for a periodic external force, the validity of LRT depends not just on the 
amplitude of the driving term but also crucially on its frequency. 

Let us consider a system characterized by a single degree of freedom, x, whose time evo- 
lution is governed by the nonlinear Langcvin equation (in dimensionless form), 

x(t) =x(t)-x 3 (t) + A cos Q,t + r){t), (1) 

where A cos fit represents an external signal and rj(t) is a Gaussian white noise with zero 
average and (rj(t)rj(s)) = 2D6(t — s). The corresponding linear Fokker-Planck equation (FPE) 
for the probability density P(x,t) reads 



f = l{(- X + X >-A C0S nt)P} + D d 1 £. (2) 
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The unperturbed system has an equilibrium distribution of the form 



D 

where N is a normalization constant and Uq(x) is the unperturbed potential 



x 4 



Uo(x)^- Y + -. (4) 

This potential has two minima located at x rn = ±1 and a maximum at Xm = 0, with 
a barrier height of 0.25. The potential Uo(x) — Axcosilt loses its bistable character for 
A > A T = 

The analysis of the dynamics is simplified by making use of two important theorems: the H- 
theorem, which ensures the existence of a unique long time distribution function Poc(x, t) ||,[| 
and the Floquet theorem, which guarantees that P oc (a;,t) is periodic in time with the same 
period as the external force 0]. For the system at hand, the symmetry of Uq(x) implies the 
following properties for the long time unique solution of the FPE: P oc (— x, t; — A) — P 00 {x 1 1; A) 
and Poo(— x,t; A) — P OQ (x,t + T/2; A), where T — 2tt/H and we have indicated explicitely the 
dependence of P^ on A. Using the Fourier expansion, 



oo 



P OQ (x,t;A) = J2 H m (x;A)e mnt (5) 

m— — oo 



the first property leads to H m (x;A) — H m (—x; —A), while the second one implies that 
H m (x; A) = (-l) m H m (-x; A). From both of them, we obtain H m (x; -A) = (-l) m H m (x; A). 
It then follows immediately that the odd moments of the distribution, (x n (t)} 00 , n = 1, 3, . . . 
can be written as Fourier series containing only odd harmonics as the even harmonics vanish 
due to the symmetries above. Analogously, even moments {x p {t))oo, p = 0, 2, . . . contain just 
even harmonics in their Fourier series expansions ||. Inserting the Fourier expansion, eq. (||), 
into the FPE, an infinite set of equations for the coefficients H m (x\ A) is obtained. Inspection 
of the set indicates that if H m (x,A) is expanded in powers of A, it can not contain powers 
smaller than j4J m L From the above general considerations, we have, in particular, for the first 
two moments, 

oo 

(^(t))^ - ^ M n {A)e innt = 2 \M n (A)\cos(nnt-cl) n ) 

n odd n>0,odd 
00 

= 2 |M n (A)| (cos(f) n cos nQt + sin^ n sinn^t) (6) 

n>Q,odd 

with M n (A) = ci 0) Al"l + c { n ] A\ n + 2 \ + ..., and 

{x(tf) x = L p (Ay pQt (7) 

p even 

with L P {A) = bi 0) AM + b { p 2) A\p+ 2 \ + .... 

The exact analytical expression for P oc (x,t) is unknown. LRT amounts to write 

Poo(M) = P eq {x) + ApW(x,t), (8) 
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with P%(x,t) obtained from a first order perturbation analysis of the FPE; see 
Appendix of [0] for details F) 



and the 



AP[ l \x,t) 



A 



n=l A « 



n 2 



[A n cos Qt + f2 sin Qt] d n (p n (x) , 



(9) 



where <p n (x) are the right eigenstates of the unperturbed FP operator and A„ the correspond- 
ing eigenvalues. The coefficients d n are ((p n \d/ dx\(po). It follows from eqs. (||) and ([)]) that 
the average value (x(t))^ iT is given by 



(x(t)) 



LRT 

CO 



oi cos(Qt - ^f^) 



The explicit calculation of the amplitude, ai, and phase lag, 



lLRT 



(10) 

, requires the knowledge of 
the spectrum of the unperturbed system. For the bistable system at hand, no exact analytical 
expressions for the eigenfunctions and eigenvalues exist, although useful approximate expres- 
sions are known @,||]. Alternatively, the amplitude and phase lag |||9|,[l(]] can be obtained 
from the response function. Using the two- mode approximation of Jung and Hanggi ||, we 
write 

(x(t))^ T = b x cos(nt - ft) + b 2 cos(fif - ft) (11) 

where the first term on the right hand side is due to the interwell hops, while the second 
one describes the influence of intrawell dynamics. It is convenient to cast the expression for 
{x(t))^ iT as in eq. ([!(]), and within the two-mode approximation, we get for the amplitude 



ai 



.4 
D 



9l\ 



x( + n 2 



1 2 

Y 2 + n 2 



2gig 2 Xia{\ia + tt 2 ) 

(\ 2 + n 2 )(a 2 + n 2 ) 



(12) 



while the phase lag of the response with respect to the input signal, < (f>i RT < 7r/2, is given 

by 

piAiO _r_ g^otO. 



arctan 



gi^ 



920S 

a 2 +n 2 



In the above formulas, Ai is given by (T^] 



Ai«^(l-§D) exp(-l/4X>), 
and a = 2. The weights, g\ and g 2 can be obtained from the expressions 



.92 



Ai(x 2 ) £ 



Ai 



Ai — a 



(13) 



(14) 



(15) 

.91 = (X 2 )eq - 92 (16) 

To leading order in D, we can replace Ai by \k = v / 2/'?'"exp(— 1/4D), g\ w 1 and 52 ~ -D/a. 
This is the limit considered in p2| . 

Linear response theory leads to the following predictions: the first moment (xft)}^ should 
contain a single harmonics with the frequency of the driving force, the output amplitude should 
behave linearly with A. Certainly, for finite values of D, if the amplitude of the driving force 



( 1 )Note that the plus signs in (A. 23) of Rcf. Jsj should read minus. This in turn yields a minus sign on the 
right hand side in (A. 28). 
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Fig. 1 - Amplitudes of the Fourier components of (a;(t))oc for noise strength D = 0.1 and input 
amplitude A — 0.2 and frequencies Q = 10 -4 (upper panel) and Q — 10 _1 (lower panel). 



is infinitcsimally small, the expansion procedure in A is valid and LRT applies. The point 
that we want to address here is that for finite small amplitudes, A < At, the value of fl has 
to be taken into account when applying LRT. The upper limit for the values of A for which 
LRT remains valid |( 2 )| depends as well on the driving frequency. 

The adiabatic approximation gives a description of the dynamics when f2 is small compared 
to any other characteristic frequency of the system. In this approach pj], the probability 
density is assumed to be given by 

^ / , T , % / Un(x) — Axcos(ilt)\ , . 

P ad (x,t)=N(t)exp( 2LJ _ L^J, (17) 

where N(t) is the normalization constant. An analysis of the corrections to the bare adiabatic 
approximation has been recently presented by Talkner . 

Even in the absence of driving, no exact explicit time dependent analytical solution of 
the FPE for the model system at hand is known. We have resorted to numerical solutions 
of eq. (||) . We follow a technique based on the use of the split propagator method of Feit et 
al. Jl5) and detailed in |ll| . From the numerical solution of the FPE we can easily obtain the 
time dependence of (x(t))oo- As this is a periodic function of time, its Fourier components 
can be obtained by numerical quadrature. 

In Fig. [l], we show the amplitudes of the relevant Fourier components of the output signal 
for D — 0.1, A — 0.2 and two very different driving frequencies, £1 = 10 _1 and = 10~ 4 . 
In this figure, as well as in the subsequent ones, we have taken D = 0.1. This is a typical 

( 2 ) In the linear response regime, the dimcnsionless ratio A/D is assumed to obey A < D. In the opposite 
singular limit, A>D, the dynamics assumes universal weak noise spectral properties |6||l3| . 
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Fig. 2 - Time evolution of (x(t)) for D — 0.1, A — 0.04 and Q. = 10 -4 as obtained from the numerical 
solution of the FPE (solid line), the adiabatic approximation (dotted line), the two- mode LRT (dashed 
line) and the two- mode LRT to leading order in D (dot-dashed line). 



value and it is adequate for the validity of the two-mode approximation leading to eqs. ( |l2| , 
|j~3| ). On the horizontal axis we indicate the order of the harmonics. It is clear that even for 
this driving amplitude, relatively large in relation to its threshold value, the response of the 
system at the larger frequency contains essentially the first harmonics. On the other hand, 
for the small driving frequency, higher order harmonics are generated. This is an indication 
of the failure of LRT to describe the dynamics at these low frequencies, while LRT might still 
be a good description for higher frequencies. 

In Fig. |5[ we depict the time evolution of (x(t)) obtained from the numerical solution of 
the FPE for D = 0.1, A = 0.04 and f2 = 10~ 4 . We also show the behaviors obtained using 
the adiabatic ansatz, eq. (|l7|), LRT within the two-mode approximation, eqs. |l4|, [Ic]) 
and LRT to leading order in D. The input signal is largely amplified at this small frequency. 
The adiabatic result deviates slightly from the numerical one at the peaks. The deviations 
from the numerical results are larger with the LRT description. Nonetheless, the two-mode 
LRT and the adiabatic approximations yield an acceptable description of the dynamics. This 
is expected within the linear response regime, where icf. 

In Fig. [| we show the behavior for fl — 10 _1 . It is clear that the adiabatic approach yields 
a signal with a very large amplitude and a large phase shift compared with the numerics. The 
two-mode LRT still yields a very acceptable behavior. The same qualitative features are 
observed in Fig. [| where £1 = 1. For this large frequency, the deviations of the two-mode LRT 
from the numerical result are very small and they can not be noticed in the plot. In Fig. ||, 
we show the full time evolution including the short transient. In Fig. ^, as we consider a large 
frequency value, we only show a few oscillations in the asymptotic regime, so that the details 
of a cycle can be distinguished. These last three figures show that for very small A, LRT gives 
a satisfactory description of the system response, with deviations from the numerics more 
pronounced as the external frequency assumes smaller values. 

To test the validity of LRT as the input amplitude is increased, we have carried out an 
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Fig. 3 - The same as in Fig. | but with Q = 10" 1 . 



extensive numerical analysis of the system response to input signals of increasing amplitudes 
and different frequencies. We evaluate the relative error e amp i = \A out — a\\/A out , between 
the output amplitude, A out , provided by the numerics and the one obtained within LRT 
with the two-mode approximation , a\ in eq. ([l2|), as a function of the input amplitude A. 
Our findings are shown in Fig. |[ The upper panel shows the dependence of e amp i on A for 
several frequencies, when the LRT is evaluated to leading order in D, while in the lower panel, 
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Fig. 4 - The same as in Fig. |2| but with Q = 1.0. Notice that due to the large value of the driving 
frequency, we only plot a few cycles of the output in the asymptotic regime. 
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Fig. 5 - Plots of the relative error of the output amplitude, e a m P i = \A out — ai\/A ut, vs. input 
amplitude A and several values of the driving frequency. In the upper panel, a\ is evaluated using 
the LRT two-mode expressions to leading order in D, while the full two-mode formulas are used in 
the lower panel; see eqs. ( |l2| , [l4[ |l5| , [l6| ) in the main text. The noise strength is D = 0.1 and the 
frequencies are: Q — 1.0 (circles), Q, = 10 _1 (plus signs), Q = 10 -3 (x) and Q = 10 -4 (triangles). 



the full expressions, eqs. ( fl2[ |TJ, [l5[ [r|) have been used. For relatively high frequencies, 
SI = 1.0 (circles), and SI = 10 _1 (plus signs), the error remains small and is practically 
constant, even for input amplitudes which are rather large compared to its threshold value. 
On the other hand, for small values of SI, fi = 10~ 3 (x) and SI = 10~ 4 (triangles), the error 
increases drastically with the input amplitude. In particular, the explicit relative errors at 
D = 0.1, i.e. (ei, e-i, £3, £4), corresponding to the driving frequencies (Sli = 10~ 4 , Sl 2 = 10~ 3 , 
fi 3 = 10 -1 , ft 4 = 1.0) respectively, read for A = 0.01: (0.028, 0.028, 0.056, 0.063); for A = 0.1: 
(0.249, 0.249, 0.072, 0.0659); and for A = 0.2: (2.539, 0.797, 0.090, 0.074). Thus, the output 
amplitude predicted by LRT at these small external frequencies is very much in error, even 
though, for the same external amplitudes and moderate-to-large frequencies, LRT predictions 
are still adequate. 

The average output lags behind the input signal with a phase shift between and tt/2. 
The value of the phase shift predicted by LRT, 4>i RT given by eqs. (|l^, [l3]), is independent 
of the driving amplitude, but depends on D and CI. It starts at for very small frequencies, 
then reaches a local maximum, and tends to its limiting value ir/2 for very large frequencies. 
Its behavior for the small-to-moderate frequencies considered here (SI < 1) is depicted with 
the solid line in Fig. || for D = 0.1. On the other hand, the phase lag of the numerical result, 

depends on D, A and SI. The \1/ values plotted have been calculated from the difference 
between the instant of times within a period, at which the driving signal and the periodic 
output, (a;^))^, cross signs, i.e., the corresponding phase delay in crossing zero. In Fig. @, we 
plot the values of ^ for several values of the driving amplitude and frequency. For very small 
frequencies, the output is almost in phase with the input for all the amplitudes considered. 
As the frequency increases, deviations between the numerical predictions and LRT results are 
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Fig. 6 - Plot of the phase lag between the average output and the driving force versus the angular 
frequency Q. With the solid line we depict 4>i m > evaluated using the LRT two-mode expressions to 
leading order in D (see eq. in the main text). The symbols denote the numerically determined 
values of the phase lag ^ (see text), for A = 0.01 (circles), A = 0.05 (triangles), A — 0.1 (diamonds), 
and A = 0.2 (x). The noise strength is set at D = 0.1. 



manifested, being larger for larger driving amplitudes. 

In conclusion, our analysis clearly indicates the influence of the driving frequency Q on the 
validity of the LRT predictions for the amplitude, phase and number of higher harmonics of 
the response of the system to subthreshold input signals. As the driving frequency assumes 
sufficiently small values, the output amplitude significantly deviates from its linear behavior 
predicted by the two-mode approximation LRT, even though the driving amplitude might still 
be quite small in order to preserve the bistable character of the unperturbed potential. Even 
for subthreshold inputs, higher order harmonics might contribute to the system response for 
small driving frequencies, contrary to the predictions of LRT. Although the global behavior of 
the phase lag indicated by LRT is qualitatively correct, as expected, its quantitative predictions 
are not reliable as the input amplitude increases. 
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